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ABSTRACT 


Previous studies have shown that terminal velocity raindrops striking a smooth water 





surface create oscillating bubbles that radiate significant underwater sound energy. Those 
studies identified two diameter ranges that produce bubbles: small drops (0.8-1.1 mm 
diameter) which produce bubbles by one mechanism and large drops (2.2-4.6 mm 
diameter) which create bubbles by a different mechanism. Effects of oblique incidence 
have been studied only for small drops. Average energy spectra were calculated for a 


range of raindrop sizes striking a smooth water surface. 


This work deals with the real life situation of large raindrops of a size often present in 
heavy rainfall (4.6 mm diameter) striking a sloped water surface. Terminal velocity is used 
to simulate natural rainfall, and the sloped surface is used to simulate the surface gravity 
waves of a natural sea. The effects of a sloped water surface on the frequency spectra and 
energy for 4.6 mm raindrops are estimated. 

By comparing energy spectra generated by single drops in an anechoic laboratory tank 
to underwater sound spectra measured at sea, it will be possible to estimate heavy rainfall 


rate by means of remote underwater listening devices. 
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I. INTRODUCTION 


One of the goals of the Naval Postgraduate School (NPS) Raindrop Lab is to 
understand the physics of underwater sound production by natural rainfall. Previous 
studies have shown that significant underwater sound is generated by the impact of 
raindrops upon a water surface and even more so by oscillating bubbles, when formed 
[Franz, 1959]. Once the sound pressure is measured, the energy spectrum of a single drop 
can be calculated. Figure 1.1 shows a typical time signal of an impact of a large diameter 
(4.6 mm ) drop which produced two bubbles . The corresponding energy spectrum is also 
shown. 

By cataloging the sound produced by each raindrop size, it is possible to predict the 
sound produced by natural rainfall [Nystuen e¢ a/ 1992]. This, in turn, will allow rainfall 
to be monitored using passive listening devices. 

Raindrops can be categorized by size and acoustic signal [Medwin ef al, 1992]. Table 
1.1 summarizes the common drop size categories and their corresponding sources of 


sound. 


TABLE 1.1. Raindrop Sizes and their Sources of Sound. 


Raindrop Diameter (mm) Sources of Underwater Sound 
i 


| Small (Type J) Impact and Bubbles 
Mid-Size 11-22 
2.2 and larger 


Large (Type II) Impact and Bubbles 
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Figure 1.1. Typical graph of a time series of on-axis acoustic pressure at 1 m for an impact and two 
resulting bubbles (upper plot) and the corresponding energy spectrum (lower plot). 
{Sampling frequency 50 kHz.; 4.6 mm drop}. 


\] 


| 
| 


The most important size categories are small and large drops as these two sizes 
consistently produce oscillating bubbles which are the dominant sound source. (The Type 
I and Type II bubble formation mechanisms are described in Chapter III.) 

All experiments at NPS have been conducted with raindrops at terminal velocity. 
Previously, however; the large drops were impacting at normal incidence upon smooth 
water surfaces. A normal angle of incidence does not describe a natural sea state. 
Although these experiments were outstanding for determining the basic physics of 
underwater sound production, a more realistic scenario is required. 

Kurgan [1989] conducted experiments for small drops at terminal velocity and at 
various angles of incidence. The influence of angle of incidence on small drops was very 
significant. His non-normal incidence angles were generated by placing a fan above the 
point of impact and forcing the small drops to strike the smooth surface at an oblique 
angle. Logistically large raindrops can not be significantly deflected using a fan. 

The purpose of this thesis is to apply the theory and knowledge gained from past 
experiments with smooth surfaces to large raindrops impacting upon a sloped water 
surface. This simulates ocean surface gravity waves. Large drops, with a diameter of 4.6 
mm at terminal velocity, were used throughout the experiment. The results will be 


compared to those for smooth surface experiments. 


Il. LABORATORY FACILITIES 


A. EXPERIMENTAL SETUP 
1. Drop Shaft 
The laboratory facilities available for this experiment are unique to the Naval 


Postgraduate School. A utility shaft with a height of 26 meters (and cross section of 3 m 





x 3 m) empties into a room containing a 1.5 m diameter, redwood lined anechoic 
cylindrical tank . [Figure 2.1] Installed in the shaft is a 9" diameter plastic tube that 
reduces the amount of interference from air drafts. The 26 m height allows drops of all 


sizes normally found in natural rainfall to reach terminal velocity. 


Intravenous Dropper 


Top View 


(for terminal 
velocity) 
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Figure 2.1. Diagram of Laboratory Setup. 
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2. Intravenous (IV) Dropper 

A standard medical IV bag with surgical tubing is used to generate large raindrops. 
Attached to the tubing is a tip calibrated to produce 50 uL drops (4.6 mm diameter). The 
accuracy of the tip was verified by measuring the volume of 100 drops six times. Each 
repetition was + 5 % of the volume expected for 100 drops. 

3. Anti-Wind Tube 

New to the lab 1s the anti-wind tube. It is constructed of thirty-three needlepoint 
hoops ( nine inch diameter) and two SO ft. lengths of typical plastic liner for gardens. The 
two lengths are overlapped by 13 feet. The hoops are equally spaced the entire length of 
the tube. The inner ring of the hoop is rolled inside of the plastic. The outer ring is made 
tight around the inner ring. This forms the plastic into a 9" x 87' cylinder. 

The addition of the tube increases the accuracy of the impact of the raindrop in 
relation to the location of the hydrophone. Previous work measured impacts with 
estimates of the random horizontal distances from the hydrophone (within 20 cm). 
Presently, we are consistently able to impact the surface of the water within 4 cm of the 
epicenter of the underwater hydrophone. The increased accuracy reduces the amount of 
range and angle correction required. 

4. Anechoic Tank 

In the room where the anti-wind tube ends is a cylindrical redwood tank with a 
diameter and a height of 1.5 m. The tank is made anechoic with a lining of redwood 
wedges. The tank contains filtered salt water acquired from the Monterey Bay Aquarium. 
Because data sets were taken on different days, the salinity was measured for each data 
set. 


The tank houses the wave generator frame and an LC-10 hydrophone. 


5. Wave Generator Frame and Motor 

Also new to the lab is the wave generator frame and motor. It is inserted into 
the tank and resembles a box kite. It is constructed of a redwood frame that supports 
walls made of Mylar. Three sides are stationary while the fourth Mylar wall is hinged at 
the bottom to allow for a paddle motion. The paddle is attached to a motor mounted level 
with the top of the tank. The dimensions of the frame are 26" x 26" x 41" (Lx Wx H). 

The wave generator remains in the tank for experiments involving both smooth and 
rough surfaces. Mylar was used because of its anechoic properties. It allows sound to 
pass but reflects surface water waves. This maintains the anechoic integrity of the tank 
and still allows for the creation of a roughened surface. 

The motor's shaft generates circular motion. This motion is translated to "linear" 


motion to drive the paddle via a mechanical eccentric coupling. 


6. Video Camera 

A Sony 8 mm video camera placed level with the water surface is used to film the 
sloped surface experiment. Extracted from the videotape are the slope of the surface and 
the depth of the hydrophone at the time of impact. The camera speed is 30 frames per 
second. 

To assist with the determination of the slope, a grid of vertical strings is suspended 
in the tank . The grid is filmed in the background while the drops impact the surface. The 
spacing between vertical lines is 2 cm. A reference line is marked on the grid before the 


surface is roughened. This allows for the measurement of the depth of the hydrophone. 


7. Hydrophone 
An LC-10 hydrophone [Figure 2.2] from Celesco Transducer Products, Inc. is 
used to measure the sound pressure of the impacts and oscillating bubbles. It is positioned 
ata depth of 6 cm for both types of experiments. For the sloped surface experiment a 
video camera is used to determine the instantaneous depth of the LC-10. The hydrophone 
IS Suspended by three supports of fishing line. The supports are positioned at intervals of 


120° to minimize the amount of movement caused by a passing wave. 


Bonded Neoprene 


Center of 
Acoustuc Field 





Figure 2.2. LC-10 Hydrophone (All dimensions are in inches). 
8. Amplifiers and Filters 


The signal from the hydrophone is connected to an Ithaco 1201 Pre-amp set at a 
gain of 100. The signal is then connected to a Krohn-Hite 3202R band pass filter. The 
frequencies passed are between |] kHz and 30 kHz in both pieces of equipment. 

9, ComputerScope® 

The signal is then connected to an IBM clone 286 computer with an analog to 
digital converter (A/D). The digital data acquisition card and its software, 
ComputerScope (sold by RC Electronics), is used for all data acquisition. 

Amplitude resolution was twelve bits. The sampling frequency used was 50 kHz. 
This allowed for a time series of 320 ms duration. Previous work was sampled at 125 


kHz and 250 kHz. Their respective time series lengths were 128 ms and 64 ms. 


Ill. BACKGROUND 


A. BUBBLE FORMATION MECHANISMS 


Two distinct bubble formation mechanisms have been identified [Snyder, 1990]. The 
Type I mechanism pertains to small diameter raindrops (0.8 mm to 1.1 mm) whereas the 
Type II mechanism pertains to large diameter raindrops (> 2.2 mm). At terminal velocity 
the raindrop diameter ranges of d< 1.1 mmand 1.1 mm<d< 2.2 mm do not produce 
bubbles. 

1. Type I 

Small raindrops at terminal velocity produce Type I bubbles 100 % of the time 
when striking a smooth surface at a normal angle of incidence. As the angle of incidence 
increases to 20°, the bubble formation percentage drops to 10 % ([Kurgan, 1989]. A 
Type I bubble is formed when the base of a conical splash crater is pinched off [Longuet- 
Higgins, 1990]. The resulting bubble resonates at approximately 15 kHz. Additional 
references for the Type I mechanism are Pumphrey ef al., [1989] and Oguz and 
Prosperetti, [1990]. This mechanism can be used to explain the sound produced 
underwater by light rain [Pumphrey ef a/, 1989] and the influence of wind upon the sound 
produced underwater by light rain [Nystuen, 1977]. 

2. Type Il 

The Type IIT mechanism is important when large raindrops are present in the rain. 
Large raindrops are prevalent during heavy rain. 

Type II bubbles are not formed as consistently as Type I bubbles. At terminal 
velocity and normal incidence upon a smooth water surface, it has been shown that 4.6 


mm drops form at least one bubble 50-65 % of the time. When Type II bubbles are 


formed they resonate between 1.6 and 10 kHz [Jacobus, 1991; Ostwald, 1992] depending 
on drop diameter. For 4.6 mm drops, the dominant resonance frequencies are between 1.6 
and 2.0 kHz. 

Snyder, Jacobus, and Ostwald describe Type II bubbles. They define a primary 
bubble as the bubble that produces the largest peak to peak voltage (pressure) for that 
drop when produced by the Type II mechanism (usually 60-70 ms after impact). 
Secondary bubbles were defined as anything else. We have found that this may not be 
appropriate. It would be better to define the "primary bubble" as the resulting bubble that 
contains the most energy. This definition is independent of the bubble production 
mechanism. Recent experiments with large drops have shown that some previously 
defined "primary" bubbles may have been produced by the Type I mechanism generated by 
late-arriving aerosols resulting from the impact. The long time after impact (>100 ms) 
when a newly defined primary bubble begins to resonate leads us to believe that it is most 
likely generated by a Type I mechanism (generated by a small drop, an aerosol, striking 
the surface). Chapter IV addresses the significant effect of a sloped surface on these "late" 
bubbles. 

Figure 3.1 (frames 1-10) shows the sequence of events that lead to the formation 
of a Type II bubble. This phenomenon was first observed by Snyder [1990]. Slow 
motion and stop motion photography of 400 frames per second allowed Snyder the 


opportunity to sketch the sequence . 


Frames 1 and 2: A flattened raindrop at terminal velocity impacts a smooth 
water surface and begins to form a canopy. The splash generates 
numerous aerosols. 


Frames 3 and 4: The canopy continues to form above the splash crater. 
Frames 5, 6, and 7: Water continues to flow up the sides of the canopy and the 


convergence of water forms upward and downward moving 
turbulent jets. 


Frames 8 and 9: The downward moving jet plunges through the bottom of 
the crater. 

Frame 10: This jet contains entrained air. A buoyant force causes the 
entrained air to break off and form a resonating bubble. 
This action occurs only if the downward moving jet is 
canted. 





Figure 3.1. Bubble Formation by the Type II Mechanism. 
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B. ENERGY EQUATIONS 
1. Conversion of Hydrophone Voltage to Pressure 


To calculate energy, the hydrophone voltage signal must first be converted to a 


pressure signal. This is done using Equation 3.1, 


V 
Pd = Go eH 


where Phyd 's the pressure at the hydrophone [Pal], V is the hydrophone voltage [volts], 


G is the amplifier gain, and ML is the hydrophone sensitivity level [volts/Pa]. 


2. Corrections to Pressure Signal 


a. Correction to One meter on Axis 
i. Dipole Radiation Pattern. Kurgan [1989] showed that the bubble 
oscillation pattern was that of a dipole (cos @ term). The geometry of the problem is 


shown in Figure 3.2. 


Location of Impact 





Figure 3.2. Location of Impact with respect to the Hydrophone. 


ii. Spherical Spreading. A correction is applied to account for the 
l/r divergence. 
Therefore, all pressures measured at range r and depth z were converted to 
I meter on axis using Equation 3.2 where rg = 1 m and cos 0 = z/r. 


Phyd 
Pim axis = ~ a (3.2) 
cos 8 1 





1] 


b. Near Field Correction (NFC) 


Because the hydrophone was within a fraction of a wavelength, an error 
between the measured pressure field and the pressure corrected to 1 meter on axis can 
occur. This can be corrected in the frequency domain using Equation 3.3 [Medwin and 


Beaky, 1989], 
(3.3) 





where & is the wave number (2n/A [m~!]) and r is the range [m] as shown in Figure 3.2. 
c. Percentage Correction Factor 
The voltage time series of a single drop was recorded if it produced at least 
one bubble. Then the energy was calculated and averaged for all drops in the data set. 
However, because only data from drops producing bubbles were taken and the fact that 
raindrops do not produce Type II bubbles 100 % of the time, the average spectral density 
was obtained by multiplying the calculated energy spectrum by a percentage correction 
factor equal to the percentage of drops which produced bubbles (determined separately). 
Jacobus and Ostwald extracted the impact and resonating bubbles from each time series. 
They then calculated the average energies due to impacts only, and then repeated the 
process for bubbles only. Jacobus showed that the impact energies are significantly less 
than bubble energies (approximately four orders of magnitude). To obtain the "average" 
spectral energy due to bubbles, the percentage correction factor was applied to the 
calculated bubble spectrum. 
This work does not separate the impacts and bubbles from the time series 
and applies the percentage correction factor to the energies calculated from the entire 


series. The error in this approach is considered insignificant. 


eZ 


| 


d. Summary 
Equation 3.4 shows all correction factors (except the percentage correction 


factor) considered at the same time for the pressure series. 


r 


lar oinaaaer field GME yp cosd NFC |V(f), G4) 
O 


V(f) is the Fourier transform of the voltage-time series. Equation 3.4 combines equations 

3.1, 3.2, and 3.3. The percentage correction factor is later applied to the ensemble energy 
density, rather than to each individual drop saved as data. 
3. Spectral Analysis and Total Energy 

A quantity proportional to the energy density spectrum can be calculated using 


Equation 3.5 [Ostwald, 1992]. 
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N is the number of points in the FFT and the frequency resolution is given by df=1/(N At). 
To convert the units to Joules/Hz, Equation 3.6 was used [Ostwald, 1992]. This 
is the energy spectral level as a function of frequency, and a subscript is added to avoid 


confusion with Equation 3.5. 
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This is the evaluation of the energy density as shown in the bottom half of Figure 1.1. 


BS 


The total energy is calculated using Equation 3.7. The subscript "s" is used to 
designate a spectral calculation and to avoid confusion with the energy (ES) which 1s 
calculated in the next section by using the measured theoretical damping constant (5) and 


the peak pressure in the time domain. 


on Tmax 
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4. Energy Calculation using the Theoretical Damping 
Constant and Peak Pressure 


A combination of equations derived by Kurgan [1989] and Scofield [1992] offers 
another method to calculate the energy of an individual bubble at a given frequency from 
the peak axial pressure and the damping constant. The combination is derived assuming 
dipole radiation and the theoretical damping constant for a given bubble size. Previously, 
Kurgan empirically calcualted the energy and did not consider the damping constant. 


Scofield was interested only in the damping constant and did not consider the dipole 





nature of a bubble. From Figure 6.3.1 of Clay and Medwin [1977], one can show that for 
the frequencies of interest for this work, the damping constant (6) can be approximated by 
5 = .0025 f (1/3) where the frequency is in Hz. The energy density on axis from an 


acoustic dipole [Scofield, 1992] can be written as 


2 
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where the pressure 1s the PEAK pressure of the bubble time series. 


The energy can be estimated by integrating Equation 3.8 (using a dipole radiation pattern) 
(Kurgan, 1989]. 


Es= \J |B cos’ @ dA (3.9) 


The elemental area is dA = (27 r sin @ )(r d@). Evaluation of Equation 3.9 yields the 


energy equation using the measured theoretical damping constant and peak pressure: 
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IV. RESULTS 


Smooth surface data was taken as well as sloped (roughened) surface data. It was 
observed that a sloped surface has a slightly higher percentage of at least one bubble being 
produced (70-73 %) than a smooth surface (50-65 %) for 4.6 mm diameter drops. The 
roughened surface data was analyzed as one data set and then again as three separate data 
sets broken into slope categories. A data set consists of a number of independent drops, 
falling at terminal velocity, and striking either a smooth or a sloped surface. A drop time 
series was saved for analysis if it produced at least one bubble. Because our time series 
extended further than ever before (320 ms), many more bubbles occurring later in time 
were detected. Jacobus’ work only extended to 64 ms and Ostwald's time series ended at 
128 ms. The duration of the series was governed by the available computer memory and 
the sampling rate of the data acquisition software (ComputerScope®). The sampling rates 
for Jacobus and Ostwald, respectively, were 250 kHz and 125 kHz. This work uses a 
sampling rate of 50 kHz. This rate is fast enough to avoid ambiguity in the signals with 
frequencies of interest and allows for a longer time series. The following presentation of 
our data will clearly show what we call the effects of a sloped water surface on the 
underwater sound radiation caused by large raindrops. For this thesis only the largest 
drop size most likely to occur in a natural, heavy rainstorm was used (4.6 mm diameter). 

All energies were calculated using MATLAB, a matrix manipulation program by 
The MathWorks, Inc. Appendix A contains the programs used to gather, process, and 
plot all data. The spectral analysis method uses Equation 3.7. The damping constant and 


peak pressure (temporal) method uses Equation 3.10. 
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A. ENERGY using the THEORETICAL DAMPING 
CONSTANT and PEAK PRESSURE 


The energy contained in dominant (primary) bubbles, as defined in earlier work 
[Jacobus, 1991; Ostwald, 1992], was calculated for the smooth and roughened surface data 
by using Equation 3.10. Previously, the energy of only the dominant bubble of each 
series was Calculated because Equation 3.4 requires knowledge of the angle @ and a range 
r from the hydrophone. Constrained by the initial definition of "primary bubble", only the 
range of the initial Type II bubble is known. Multiple hydrophones will be needed to 
accurately measure the range to all bubbles created. However, the use of the damping 
constant and the peak pressure to evaluate the energy for each bubble allows for an 
estimate of the energy in the entire series. 

1. Energy at Higher Frequencies 
Figure 4.1 (top) shows the same general distribution as reported by Ostwald. 
Figure 4.1 (bottom) shows that for a sloped surface at the time of impact, there exists 
more energy at higher frequencies. Higher frequency bubbles are produced by smaller 
raindrops [Kurgan, 1989; Jacobus, 1991; Ostwald, 1992]. Because they are higher 
frequency, one would assume that they are formed by the Type I mechanism (or late 
aerosol impacts). There is no definitive proof that this is the case. The roughened 
surface could possibly produce Type II bubbles with frequencies higher than previously 
observed. 
2. Time Gap 
As previously stated, earlier work used a shorter time series. When the 
duration of the series is extended, a noticeable time gap was observed for the smooth 
surface data. This gap occurs between 110 and 140 ms as suggested in Figure 4.2 (top). 


This time gap is not present in the roughened surface data [Figure 4.2, bottom]. 
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Figure 4.1. Comparison of Energies of Dominant Bubbles for Smooth (top) and Roughened 
(bottom) surfaces. 
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Figure 4.2. Comparison of Dominant Bubble Energy vs Time after Impact for 


Smooth (top) and Roughened (bottom) surfaces. 


We 


When the frequencies of all bubbles that were present in the time series (primary 
and secondary) are plotted against the time after impact, the time gap becomes more 
prevalent in the smooth surface data [Figure 4.3, top]. The end of the gap (approximately 
140 ms after impact of the large drop) corresponds to the time delay for an aerosol 
trajectory with a maximum height of 2.5 cm. Bubbles occurring after 140 ms would not 
have been considered as dominant bubbles in previous studies. Because of the timing as 


shown in Figure 3.1, they can not be Type II bubbles. 


Aerosols generated by impacting a smooth, stationary surface tend to be 
"launched" almost vertically and have been seen reaching extreme heights of 40 cm above 
the surface. This corresponds to a time of flight (time after impact) of 571 ms. Visually, 
the average height was between 12 and 15 cm. These heights correspond to times after 
impact of 312 and 350 ms, respectively. This suggests that the duration of the data 
collection should be lengthened further, although there would be a corresponding 
reduction in sampling rate. 

For the roughened surface data in Figures 4.2 and 4.3 (bottom) the time gap 1s 
not present. This is consistent with a reduced time of flight caused by a lower trajectory 
of the aerosols. The rough surface aerosols are also launched with a much greater 
horizontal velocity component than the smooth surface aerosols. It was quite obvious that 
the horizontal velocity of the surface wave, as well as the wave slope, greatly affect the 


velocity components and the launch angles of these aerosols. 
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Figure 4.3 Comparison of Frequencies of ALL Bubbles vs Time after Impact for 
Smooth (top) and Rough (bottom) surfaces. 
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B. SPECTRAL ANALYSIS 
1. Smooth vs Roughened Surfaces 

Figure 4.4 compares the average spectral energy density of the smooth surface 
data (top half) to the average spectral energy density of the roughened surface data 
(bottom half). The instantaneous surface slope at the time of impact was not measured 
because we were looking for a general effect due to a randomly sloped surface. It can be 
seen from Figure 4.4 that the average total energies are comparable. This was not 
expected because more high frequency bubbles were observed in the roughened surface 
data compared to the smooth surface data. Typically, high frequency bubbles radiate less 
energy than low frequency bubbles; therefore their influence on the total energy is not as 
significant as lower frequency bubbles. 

The plot of the average energy density for the roughened surface in Figure 4.4 
(bottom) also demonstrates an apparent shift of the peak of the spectrum to higher 
frequencies. This was expected because of the presence of higher frequency bubbles 


occurring late in the time series. 
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Figure 4.4. Comparison of Average Energy Densities for Smooth (top) and Roughened 
(bottom) surfaces. 
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2. Roughened Surface Categorized by Slope 

The use of a video camera allowed us to determine the slope of the surface at the 
time of impact. The slopes were grouped into three categories. The mean slopes were 0, 
6, and 11° from the horizontal (t 3°). A 0° slope for a roughened surface is different 
from a smooth surface because of the horizontal velocity of the surface waves. The wave 
velocity was seen to be imparted to the raindrop at impact and affected the nature of the 
aerosols generated by the impact. 

Figure 4.5 shows the three slope categories. It can be seen that the average total 
energies are again comparable, independent of the slope of the surface at the time of 
impact. Forty-eight drops were used for this analysis. The energies are all within one 
standard deviation of each other. Better statistics are needed before stronger statements 
can be made. 

As with the randomly roughened surface of Figure 4.4, the peaks of the spectra of 


Figure 4.5 also shift to higher frequencies as the slope of the surface increases. 
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Figure 4.5. Comparison of average energy densities of three surface slope categories. 
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V. CONCLUSIONS 


The original purpose of this thesis was to determine the effects of a sloped water 
surface at the time of impact for a 4.6 mm diameter raindrop. While several effects were 
observed, this experiment also opens the door to numerous other studies. The results of 


this study reveal that: 


@ For smooth and rough surfaces (including slope categories): average energies per 
drop are of the same order of magnitude [Figure 4.4 and Figure 4.5]. 

@ The peak of the energy spectrum for the sloped surface apparently shifts to higher 
frequencies [Figure 4.4, bottom and Figure 4.5]. 

6 The "time gap" discovered for smooth surface data is absent for the sloped surface 
data. We now understand that the time gap is absent in the sloped surface data 
because the characteristics of sloped surface aerosols are greatly affected by the 


slope of the surface and the horizontal velocity of the surface wave. 


Only after more is known about the physics of underwater sound produced by 
raindrops of different sizes generated in a laboratory, will we be able to understand the 
physics of underwater sound produced by natural rainfall at sea. This will ultimately 
establish the ability to estimate rainfall rate by means of remote underwater listening 


devices. 
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VI. RECOMMENDATIONS 


A. DEFINITIONS 


1. Primary 
Redefine Primary (dominant) to be the most energetic bubble. This will 
remove any dependence on its origin (Type I or II). 
2. Type I Mechanism 
Only include bubbles created by pinch-off from conical cavity. 
3. Type II Mechanism 
Only include bubbles created by downward jet within the canopy of large drop 
splashes. 
B. FUTURE EXPERIMENTS 
1. More Data for this Experiment 
In order to establish respectable statistics, more data needs to be taken. More 
data will also help to enforce the results of this work. The effect of a sloped surface for 
different raindrop sizes should be studied. 
2. Multiple Raindrops 
So far only single raindrops have been used for all experiments. Design a 
multiple drop experiment to determine if there is any interaction between the drops at 
impact (linear or nonlinear). The number and volume of the drops must be controllable in 


order to better understand the resulting spectra. 


Zz 


3. Multiple Hydrophones 
Multiple hydrophones can be used to determine ranges of impacts or bubbles. 
This will determine the origin (bubble production mechanism) of all bubbles occurring 
later in time. 
4. Capillary waves 
Generate a roughened surface by creating capillary waves and then perform 
this same experiment. Compare the results of the smooth, sloped, and roughened 
(capillary) surface data. 
5. Eliminate the % Correction Factor 
Record the data series even if bubbles are not present. This will allow for a 
more accurate average total energy for a known number of drops and bubbles. 
6. Extend the Listening Time 
Use time delays of 320 ms or 640 ms to extend the time series to 640 ms or 


960 ms, respectively. This will not affect the sampling rate. 
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APPENDIX A 


The following three programs were written to analyze the single drop data collected 
in the NPS raindrop tank. The data was collected using ComputerScope© and stored in 
an ASCII file format. These files were stripped of extraneous headers and time 
information and stored as data files containing only the output voltage. (The time step is 
known and is input inside these programs). The stripped ASCII files are imported into 
MATLAB. These programs are written in the MATLAB processing language 
(MATLAB *.m files). All plotting routines can be easily changed within the program to 
present data any way the user desires. 

The first program, ANAL.M, extracts: 

the frequency of any bubble existing in the time series, 
the time after impact of each bubble, 
and the peak to peak voltage of each bubble in the time series. 
The data is broken up into "dominant" and "all". It then saves the extracted data in 


matrices used in succeeding programs. 


% Program ANAL.M 

Yo Program to ANALyze data taken by ComputerScope and converted to ASCII. 
‘% Calculates the frequency, time after impact, and Vpp of a particular bubble. 
Yo Two mouse "ginputs" are required by the user. This is used for the freq. 

‘% It also creates and saves (user chooses name of file) the data in a *.mat file 
Yo for further manipulation in *.m files like delta.m. 

ele 

clear 

clg 

disp(‘In what directory will I find the data ‘') 

direct = input(‘{i.e. full path [in single quotes]} ? ....'); 

disp(’ ') 

disp(’ ’) 

start = input(‘Input number of starting file...'); 

stop = input(‘Input number of ending file...'); 


Zo 


maxnum = stop - start + 1; 


ftV=| |: % empty matrix to allow the building of the true data matrix 
ftVdom={ J; % empty matrix to allow the building of the true data matrix 
cle 


disp(‘With what letter does the data file name start ') 
a=input(‘{i.e. r*.dat, s*.dat, etc} [in single quotes]...?'); 
fori = l:maxnum; 

clg 

cle 

eval({'load ',direct,'\'",a,int2str(start + 1 - 1), '.dat')) 

eval({'v = ',a,int2str(start + 1 - 1),';')) 

eval({‘clear ',a,int2str(start + 1 - 1),";'}) 

delt=20e-6; 

vV =v - mean(v); 

plot(v) 

xlabel(‘n') 

ylabel(‘{volts]') 

title({"‘Drop # 'num2str(i),' for this data set.']) 

grid 

pause 

cle 

disp('Entering keyboard mode so that the user can determine the # of’) 

disp(‘bubbles for this time series (in case of uncertainty).’) 

disp(' ') 

disp(‘[{CTRL-z] to exit keyboard mode’) 

disp(’ ’) 

disp(’ ') 

keyboard 

ele 

bub=input('How many bubbles in this series?...'); 


for b=1:bub; 
cle 
low=input(‘Input lower data point (n) for plotting...'); 
up=input(‘Input upper data point (n) for plotting...'); 
subplot(211),plot(v) 
xlabel(‘n') 
ylabel(‘{volts]') 
title({‘Drop # ',;num2str(i),' for this data set']) 
gnid 
subplot(212),plot(v(low:up)) 
xlabel(‘n') 
ylabel(‘[volts]") 
title({"‘Bubble # ';num2str(bub),' for this drop']) 
grid 
pause 
disp(’ ’) 
check=input(‘Are these the correct limits?...1-Yes 2-No...’); 
while check == 2 
cle 
disp(' ’) 
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disp(‘Last limits were... ‘) 
disp(’') 
limits=[low up] 
disp(’ ’) 
low=input(‘Input lower data point for plotting...'); 
up=input(‘Input upper data point for plotting...'); 
clg 
subplot(211),plot(v) 
xlabel('n') 
ylabel(‘{[volts}') 
title(["Drop # ';num2str(i),' for this data set']) 
grid 
subplot(212),plot(v(low:up)) 
xlabel(‘n’) 
ylabel(‘[volts]') 
title({"Bubble # ',num2str(bub),' for this drop']) 
grid 
pause 
clear check 
check=input(‘Are these the correct limits?...1-Yes 2-No...’); 
end % end while 
cle 
[mx,c]=max(v(low:up)); 
mn=min(v(low:up)); 
[xt yv]=ginput(2); 
freq=(abs(1/((xt(1)-xt(2))*delt)))/1000 % [kHz] 
=((ctlow)*delt)* 1000 % [ms] 
Vpp=abs(mx) + abs(mn) % [volts] 
disp(’ ') 
ftVb=[freq t Vpp]; % freq, time, Vpp for each bubble 
disp(’ ') 
if bub~=1; 
dom=input('Is this the dominant bubble of the series?...1-Yes 2-No...'); 
else 
dom=1; 
end 
if dom==1; 
ftVdom=[ft Vdom;ft Vb]; 
end 
ftV=[ftV;ftVb]; 
clg 
plot(v) 
title({'Extra look at drop # ';num2str(i),' for this data set']) 
grid 
pause 
cle 
end % next b 
cle 
clear v 
end % next i 
mat=input('What name for the *.mat file [in single quotes] ? ....'); 
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disp(‘ ') 

disp(‘The matrix containing...freq [kHz] t [ms] Vpp [volts]...will be saved’) 
disp(‘for future uses (i.e. delta.m). The variables will be ftV and ftVdom.") 
disp(’') 

disp(‘Insert a disk into drive a: to receive the *.mat file...') 

disp(‘[ENTER] to continue’) 

pause 

la: 

eval([‘save ',mat,' ftV ftVdom']) 

Ic: 

cle 

disp(‘It is now recommended to set up and run DELTA.M') 
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The second program, DELTA.M, uses the matrices formed in ANAL.M and 
calculates the theoretical damping constant (6). Then the energy using 6 and the 1/2 of 
Vpp ( peak voltage) found in ANAL.M is calculated. Finally, many different combinations 


of data are plotted like Figures 4.1, 4.2 and 4.3. 


% Program DELTA.M - 
% to be used after running ANAL.M and before ANAL2.M 


‘“ Creates monster matrix and calculates the damping constants and 
% Energy based on these deltas (Clay & Medwin text Ch 6) 
% Plots different combinations of the data. 


format short e 

cic 

clear 

pack 

clg 

disp(‘Data consists of the matrices created in ANAL.M. It is') 
disp(‘the *.mat file containing the matrices consisting of:’) 
disp(’ ') 

disp(‘[ f (kHz) t (ms) Vpp (volts) J.') 

disp(’ ’) 

disp(’ ') 

disp(' What was the name of the *.mat file from’) 

file = input(ANAL.M {FULL PATH including filename (in single quotes)}.....?'); 
cle 

disp(’Where are the files d.dat and z.dat') 

direct2 = input('{i.e. FULL PATH (in single quotes)}.....?'); 
eval([‘load ',file]); 

eval([‘load ',direct2,'\’,'d.dat']); 

eval([‘load ',direct2,'\’,'z.dat']); 


d=d./100; % [m] 
z=z./100; % [m] 
d=d(1:length(ftVdom(:,!))); % [m] 
z=z(1:length(ftVdom(:,1))); % [m] 


% thoc = input(‘Input the Specific Acoustic Impedance (rhoc)...'); 
“ diam = input(‘Input drop diameter [mm]...'); 

% gain = input(‘Input hydrophone gain...'); 

% ML = input(‘Input hydrophone Sensitivity [V/Pa]...'); 

r0=1; % 1 meter on axis 

c = 1500; % sound speed 

rhoc = 1.54e06; % for seawater 

ML = 2.66e-5; % Hydrophone sensitivity (V/Pa) over 1-20 kHz 
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diam=4.6; 
gain=100; 


% Calculate the damping constant based on a linear regression of the 
% delta vs. freq (log-log plot) from Clay & Medwin 


del=.025 * ftVdom(:,1).%(1/3); % ftVdom(:,1) is freq column [kHz] 
“% Calculate energy based on damping constant ‘delta’ 


r=sqrt(z.“2 + d.*2): “% [m] 

k=2*pi*ftVdom(:, 1)*1000/c; % 1000 to convert to Hz 
CF=(r.42)./(gain*ML*1r0*z); 

NEC = (k.*r) 7 (sqrt((ki*1).*> 24 

pax=NFC. *CF.*ftVdom(:,3)/2; % 1 mon axis based on PEAK voltage 
E=((r.*pax).42)./(1000*6*rhoc.*ft Vdom(:,1).*del); % 1000 to convert to Hz 


ft VdEdom=[ftVdom del E ]; 


% This part of delta.m plots different combinations of ftVdEdom. 

% This matrix now consists of five columns: 

% [f(kHz) t(msec) Vpp(Volts) delta(damping constant) Energy(joules)] 
cle 

disp(‘What type of surface was used for this data’) 

surf = input('{i.e. Smooth or Roughened (in single quotes)}.....?'); 

cle 


name = input('What name do you want for the *.met file {in single quotes}..... 


cle 

disp(‘What name do you want for the *.mat file containing ') 

name2 = input(‘the [ ft Vpp delta E } matrix {in single quotes}.....? '); 
cle 

m=mean(ft VdEdom); 


% plot #1 

“ Dominant bubble Energy vs freq 

% 11=[0, 25, -12, -6]; 

%  axis(11) 

%  subplot(211),semilogy(ftVdEdom(:,1),ftVdEdom(:,5),'+') 
%  tit=["Dominant Energy vs Frequency...',surf,' Sfc']; 

%  title(tit) 

%  xilabel('freq [kHz]}') 

%  ylabel('Energy [Joules]') 

% gtext(['Avg Energy = ',num2str(m(5)),' Joules')) 


% plot #2 

“% Dominant bubble Energy vs time 

% 12=([0, 300, -12, -6}; 

% axis(12); 

%  subplot(211),semilogy(ftVdEdom(:,2),ftVdEdom(.:,5),'+') 
%  subplot(212),semilogy(ftVdEdom(:,2),ftVdEdom(.,5),'+') 
%  tit=['Dominant Energy vs Time after Impact...',surf,' Sfc'); 
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%  title(tit) 
%  Xlabel(‘time [msec]') 
%  ylabel(‘Energy [Joules]') 


% cle 

% pause 

% disp(‘Insert a disk into drive a to receive the *.met and *.mat files.’) 
% disp((ENTER to coninue’) 

“% pause 

% ‘a: 

% eval(['‘meta';name]) % plots #1 & #2 

% 'C: 

% cle 


% plot #3 
“% Dominant bubble Energy vs drop # 


%  subplot(211),semilogy(ftVdEdom(.,5),'+') 

%  tit=('Dominant Energy vs Drop number...',surf,’ Sfc']; 
%  title(tit) 

%  xlabel('Drop #') 

%  ylabel('Energy [Joules]') 


% plot #4 
‘ All bubbles...Freq vs time 


14=(0, 300, 0, 25]; % force the axis limits for direct comparisons 
axis(14) 
subplot(211),plot(ftV(:,2),fV(-, 1),'+') 
gtext('ALL bubbles’) 
%  subplot(212),plot(ftV(:,2), VC, 1),'+') 
tit=['"Frequency vs Time after Impact... ‘surf,’ Sfc']; 
title(tit) 
xlabel('t [msec]') 
ylabel(‘freq [kHz]') 


pause 

% !a: 

eval(('meta 'sname]) 

eval({'save 'sname,' ft VdEdom']) 

% 'C: 

cle 

disp(’ ') 

disp(’ ’) 

disp(‘It is now recommended that you set up and rin ANAL2.M’) 
disp(’ ’) 

disp(‘Thank you for playing DELTA..') 


Sie 


The final program, ANAL2.M, calculates and averages spectra for a series of drops 
and produces plots like those shown in Figures 4.4 and 4.5. 


% Program ANAL2.M 


“% for rough and smooth surfaces 


Y% program to ANALyze spectral level determination of raindrop data 
‘“o Converts bubble time series to power spectrum. Outputs result. 


cle 

clear 

clg 

c = 1500; % sound speed (use the correct sound speed for the salinity) 
rhoc = 1.54e06; % for seawater 

ML = 2.66e-5; % Hydrophone sensitivity (V/Pa) over 1-20 kHz 

diam=4.6; 

gain=100; 

N=16384; 

delt = 20 * 10“(-6); 

disp(’ ') 

disp(' NOTE: Drop file names must be of the form r#.dat or s#.dat.’) 

disp(' The drop files must be numbered consecutively.') 

disp(' ') 

disp(’ Ensure impact distances are in file d.dat. ') 

disp(’ and hydrophone depths are in file z.dat.') 

disp(’ D.DAT and Z.DAT MUST be in the same directory as your data files.’) 
disp(’ Distances must be in a single column, in cm’) 

disp(’ ') 


dcheck = input('Did you update d.dat (distance) and z.dat (depth) files (1 = yes, 2 = no)..'); 
if dcheck ~= 1 
disp(’ ') 
error(‘Update the files and rerun this program.') 
end 
clear dcheck; 


cle 


% diam = input(‘Input drop diameter (mm)...'); 

Y gain = input(‘Input hydrophone gain...'); 

‘“ ML = input(‘Input hydrophone Sensitivity [V/Pa]...'); 

‘ N = input(‘Input number of points for FFT...'); 

‘% delt = input(‘Input sample period (time between samples, usec)..."); 
Y% delt = delt * 10“(-6); 

pack 

cle 
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surf = input('What type of surface? (1-Smooth 2-Roughened) '); 
disp(' ') 


if surf==2 % Rough Surfaces 
cle 
disp((DATA includes z.dat, d.dat. and r*.dat files’) 
disp(’ ’) 
disp(‘In what directory will I find the DATA’) 
direct = input(‘{i.e. full path [in single quotes]} ? ...'); 
eval([‘load ‘direct ,'\d.dat;']) 
d=d./100; % [m] 
eval([‘load ‘direct ,'\z.dat;']) 
=z./100; % [m] 
cle 
name = input(‘What do you want to call the meta file [in single quotes] ? ...'); 
disp(’ ') 
disp('What is the Sfc slope for this data (i.e. "5.5 degrees")’) 
slope = input(‘[in single quotes] ?...'); 
cle 
disp(‘For Roughened Surfaces...') 
start = input(‘Input number of starting file...'); 
stop = input(‘Input number of ending file...'); 
maxnum = stop - start + 1; 
perc = input(‘Input % of drops that produce at least one bubble...'); 
perc = perc / 100; 
df = 1/(N * delt); 
Esum =zeros(1:(N/2)); % Initializing size of row vector 
f= (df.* [0:1:(N/2)-1])'; 
ke=1@ * pimted)./c; 
r0=1; % [m] 


for1 = l:maxnum; 
eval({‘load ',direct,'\r',int2str(start + 1 - 1), '.dat']) 
eval(('v =','r',int2str(start + 1 - 1),";']) 
eval({'clear ','r',int2str(start + 1 - 1),';']) 
v =v - mean(v); 
% t=0:delt:(length(v)-1)*delt; 
V = abs(fft(v,N)); 


clear v 
r(i) = sqrt(z(i* 2 + dCi). 2); % [m] 
NFC = (k.*r(i)) ./ (sqrt((k.*r(i)).* 2 + 1); 
CF=((perc*r(1)*2)/(gain*ML*z(i)*1r0)); % [Pa/v] 
Pax = CF .* NFC .* V(1:(N/2)); % [Pa] 1 mon axis 
% Using eqn 41 from thesis by Ostwald 
Ef=((2*delt)./(N*df)).* Pax(1:N/2).%2; % [Pa*2-s/Hz] 
clear Pax V 
E=((2*pi*r(i).42)/(3 *rhoc))*Ef, % [Joules/Hz] 
clear Ef 
Esum = Esum + E'; % for averaging purposes later 
clear E 
end 
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clear CF NFC start stop N delt i surf r 


Eavg = Esum/maxnum; 
etot=sum(Eavg) *df; 


SL = (10*log10(Esum) + 120)'; % dB re 1 uPa‘*2/Hz 


cle 


disp(‘Four gtexts follow: Sfc slope, Avg Energy, slope dependence, freq shift’) 
disp(‘[Enter to continue]') 
pause 
11=[0, 18, 0, 4e-13}]; 
axis(11) 


% 
Yo 


subplot(211),plot(f./1000,Eavg); 
subplot(212),plot(f./1000,Eavg); 


plot(f./1000,Eavg); 
xlabel('freq (kHz)') 
ylabel('[Joules / Hz]') 
if maxnum==1; 


title({'Avg Energy Density for ';num2str(maxnum),' drop for Roughened Sfc']) 


else 


title(['Avg Energy Density for ';num2str(maxnum),' drops for Roughened Sfc’}) 


end 

%('Drop Diameter = ',;num2str(diam),' mm']) 
gtext(['Sfc slope = ',|slope]) 

gtext(['Average Total Energy = ';num2str(etot),’ [J]']) 
gtext(‘No significant dependence on slope’) 
gtext(‘Spectrum shifts to higher frequencies’) 


% 
% 
% 
% 
% 
% 
0 


% 


S 


subplot(2 12), plot(f(2:8190)./1000,SL(2:8190)); 
xlabel(‘freq (kKHz)') 
ylabel('dB re 1 uPa%2 / Hz’) 
if maxnum==1; 

title(['Spectrum Level for ';num2str(maxnum),' drop']) 
else 

title(['Spectrum Level for ';num2str(maxnum),' drops']) 
end 


pause 


Yo 
% 


cle 

disp(‘Insert the *.met disk into A’) 
disp(‘[Enter to continue]') 

pause 

rae 


eval(['meta ';name]) 


% 


le: 
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% Average Energy density of all drops 
% Evaluate the integral eqn 42 


else % Smooth Surfaces 


cle 

disp(’DATA includes d.dat, z.dat, and s*.dat files') 

direct = input(‘In what directory will I find the data {1.e. full path [in single quotes]} ? ....'); 
eval([‘load ',direct,'\d.dat']) 

eval([‘load ',direct,'\z.dat']) 

d=d./100; % [m] 

z=z./100; % [m] 

name = input(‘What do you want to call the meta file [in single quotes]} ? ...'); 
cle 

disp(‘For Smooth Surfaces...') 

start = input(Input number of starting file...'); 

stop = input(‘Input number of ending file...'); 

maxnum = stop - start + 1; 

perc = input(‘Input % of drops that produce at least one bubble...'); 
pere — perc / 100; 

df = 1/(N * delt); 

Esum =zeros(1:(N/2)); % Initializing size of row vector 
f= (df.* [0:1:(N/2)-1))’; 

ke (2*-pi.* f)./c 

fe sart(z. 2 +d.“ 2); % [ml] 

r0=1; % [m] 


for 1= l:maxnum; 
eval([‘load ',direct,'\s',int2str(start + 1 - 1), '.dat']) 
eval(['v =','s',int2str(start + i - 1),';']) 
eval([‘clear ','s',int2str(start + i - 1),';']) 
vV=v- mean(v); 
% t=0:delt:(length(v)-1)*delt; 
V = abs(fft(v,N)); 
clear v 
1 “% Progress Pointer for the user 
NFC = (k.*r(i)) ./ (sqrt((k.*r(i)).* 2 + 1)); 
CF=((pere*r(i)*2)/(gain* ML*z(i)*10)); % [Pa/v] 
Pax = CF .* NFC .* V(1:(N/2)); % [Pa] 1 monaxis 
% Using eqn 41 from thesis by Ostwald 
Ef=((2*delt)./(N*df)).* Pax(1:N/2).%2; % [Pa*2-s/Hz] 
clear Pax V 
E=((2* pi*r(i).42)/G *rhoc))*Ef: % [Joules/Hz] 
clear Ef 
Esum = Esum + E'; % for averaging purposes later 
clear E 
end 


clear CF NFC start stop N delti surf r 


Eavg = Esum/maxnum; % Average Energy density of all drops 
etot=sum(Eavg)*df; % Evaluate the integral eqn 42 
SL = (10*logl0(Esum) + 120)’; % dB re 1 uPa%2/Hz 


a0 


disp(’ ') 
%  disp('Please insert the *.met disk into drive A [Enter to continue]') 
Y% pause 

cle 

disp(‘One gtext follows: Avg Energy’) 

disp(‘{Enter to continue]') 

pause 


11=[0, 18, 0, 4e-13]; 
axis(11) 
% subplot(211),plot(f./1000,Eavg); 
plot(f./1000,Eavg); 
xlabel(‘freq (kHz)') 
ylabel('[Joules / Hz]') 
if maxnum==1; 
title({'Avg Energy Density for ',num2str(maxnum),' drop for Smooth Sfc']) 
else 
title({/Avg Energy Density for 'snum2str(maxnum),' drops for Smooth Sfc’]) 
end 
“gtext({'Drop Diameter= ',;num2str(diam),' mm']) 
gtext(['Average Total Energy= ',num2str(etot),’ [J]']) 


% subplot(212),plot(f(2:8190)./1000,SL(2:8190)); 

% xlabel('freq (kHz)') 

% ylabel('dB re 1 uPa‘%2 / Hz’) 

% if maxnum=1; 

%  title({'Spectrum Level for 'snum2str(maxnum),' drop']) 
% else : 
o  title({'Spectrum Level for ';num2str(maxnum),' drops']) 

% end 














SS 


pause 
Yourla: 

eval(['meta ',name;]) 
% Ic: 
end 
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